#include "algo-support-reconnect.h"
#include <list>
#include <boost/math/special_functions/factorials.hpp>
#include <boost/math/special_functions/binomial.hpp>
#include <boost/foreach.hpp>
#include <boost/tuple/tuple.hpp>

#define foreach         BOOST_FOREACH
#define reverse_foreach BOOST_REVERSE_FOREACH

using namespace std;

namespace cccn
{
namespace internal
{

////////////////////////////////////
//                                //
// Definition of internal classes //
//                                //
////////////////////////////////////

MultinaryCounter::MultinaryCounter () : m_counter (0), m_weight (1,0), m_max_counter (0)
  {
  Reset ();
  }

MultinaryCounter::MultinaryCounter (vector<uint> max) : m_counter (0), m_weight ()
  {
  Reset (max);
  }

void 
MultinaryCounter::Reset ()
  {
  m_counter = 0;
  }

void 
MultinaryCounter::Reset (vector<uint> max)
  {
  Reset ();

  // initialize the m_max_counter and m_weight
  m_max_counter = 1;
  m_weight = vector<uint> (max.size (), 1);
  vector<uint>::reverse_iterator weight_it = m_weight.rbegin ();
  reverse_foreach (uint m, max)
    {
    // initialize the m_weight for each element. the weight of the
    // last element must be 1
    *weight_it = m_max_counter;
    weight_it++;
    // calcluate the maximum possible number for the m_counter
    m_max_counter *= m;
    }
  }

bool 
MultinaryCounter::operator++ (int)
  {
  if (m_counter+1 < m_max_counter)
    {
    m_counter++;
    return true;
    }
  else
    {
    return false;
    }
  }

uint 
MultinaryCounter::at (uint i) const
  {
  if (i == 0)
    {
    return m_counter / m_weight.front ();
    }
  else if (i-1 == m_weight.size ())
    {
    return m_counter % m_weight[i-2];
    }
  else
    {
    return (m_counter % m_weight[i-1]) / m_weight[i];
    }
  }

/////////////////////////////
//                         //
// Definition of functions //
//                         //
/////////////////////////////

vector<uint>
GetPermutationAt (const uint n, const uint r, uint i)
  {
  // store the resulting pattern
  vector<uint> pattern (r, 0);
  // initialize all available to be drawn
  list<uint> items;
  for (uint x=0; x<n; x++)
    {
    items.push_back (x);
    }

  for (uint x=0; x<r; x++)
    {
    const uint skip = boost::math::falling_factorial (n-x-1, r-x-1);
    // initially, pattern stores not the actual pattern, but the index for drawing an item
    uint draw_at = i / skip;
    // upper bound the parameter i
    if (draw_at >= items.size ())
      {
      draw_at = items.size ()-1;
      }
    // go to the item indexed by 'draw_at'
    list<uint>::iterator draw_it = items.begin ();
    for (uint counter=draw_at; counter>0; counter--)
      {
      draw_it++;
      }
    // draw the item and store it into the 'pattern'
    pattern[x] = *draw_it;
    // remove the item drawn from the items list
    items.erase (draw_it);
    // 
    i -= draw_at * skip;
    }

  return pattern;
  }

vector<uint>
GetCombinationAt (const uint n, const uint r, uint i)
  {
  vector<uint> pattern (r, 0);
  
  uint draw_item = 0;
  for (uint x=0; x<r; x++)
    {       
    for (uint y=draw_item+1; n-y>=r-x-1; y++)
      {
      // upper bound the i
      if (draw_item >= n-(r-x))
        {
        break;
        }
      const uint skip = boost::math::binomial_coefficient<double> (n-y, r-x-1);
      if (i >= skip)
        {
        i -= skip;
        draw_item++;
        }
      else
        {        
        break;
        }
      }
    // construct the pattern
    pattern[x] = draw_item;   
    draw_item++;
    }
  return pattern;
  } 

////////////////////////////////////////////////////
//                                                //
// Definition of ReconnectionPatternCounter class //
//                                                //
////////////////////////////////////////////////////

ReconnectionPatternCounter_odd_k::ReconnectionPatternCounter_odd_k (const Graph& g) 
: m_g (g), m_edge_pattern_index (), m_degree (0), m_v1 (0), m_v2 (0)
, m_reconnect_root (), m_node_pattern_index (0), m_reconnect_edge_index (0)
, m_maxNodePatterns (0), m_index (0)
  {
  typedef boost::graph_traits<BoostGraph_t>::adjacency_iterator adjacency_iterator;    

  const BoostGraph_t& boostg = g.GetBoostGraph ();
  const uint n = num_vertices (boostg);

  // the two new nodes can be found easily by looking for the largest and the second largest node IDs
  m_v1 = n - 1;
  m_v2 = m_v1 - 1;
  Vertex_t boost_v1 = vertex (m_v1, boostg);
  Vertex_t boost_v2 = vertex (m_v2, boostg);
  // initialize the m_degree, which is the expected degree of every node after reconnection
  m_degree = out_degree (boost_v1, boostg) * 2 - 1;

  /* initialize m_reconnect_root
   *
   * note about the structure
   * {... va_adj ...} === {... va ...} --- v1 --- v2 --- {... vb ...} === {... vb_adj ...}
   * reconnection are done by changing one of the edges 
   * va_adj---va to va_adj---v1 or va_adj---v2.
   * The same logic applies to nodes vb_adj
   */
 
  Vertex_t boost_v1p, boost_v2p;
  for (uint i=0; i<2; i++)
    {
    if (i == 0)
      {
      boost_v1p = boost_v1;
      boost_v2p = boost_v2;
      }
    else
      {
      boost_v1p = boost_v2;
      boost_v2p = boost_v1;
      }

    // for each va (or vb), obtain its neighboring nodes and keep them into the m_reconnect_root
    adjacency_iterator boost_v1p_adj_it_begin, boost_v1p_adj_it_end;
    tie (boost_v1p_adj_it_begin, boost_v1p_adj_it_end) = adjacent_vertices (boost_v1p, boostg);
    for (adjacency_iterator va_it = boost_v1p_adj_it_begin; va_it != boost_v1p_adj_it_end; va_it++)
      {    
      // skip the v1 and v2
      if (*va_it == boost_v2p)
        {
        continue;
        }
      uint va_id = get (boost::vertex_index, boostg, *va_it);
      map<uint, vector<uint> >::iterator reconnect_root_it;
      bool vaIsNew;
      boost::tie (reconnect_root_it, vaIsNew) = m_reconnect_root.insert (make_pair (va_id, vector<uint> ()));

      if (vaIsNew)
        {
        // the v_id does not exist in the m_reconnect_root vector

        // obtain all neighbors of va
        adjacency_iterator va_adj_it_begin, va_adj_it_end;
        tie (va_adj_it_begin, va_adj_it_end) = adjacent_vertices (*va_it, boostg);
        for (adjacency_iterator va_adj_it = va_adj_it_begin; va_adj_it != va_adj_it_end; va_adj_it++)
          {
          // skip the v1 and v2
          if (out_degree (*va_adj_it, boostg) != m_degree)
            {
            continue;
            }
          uint va_adj_id = get (boost::vertex_index, boostg, *va_adj_it);
          // store the va_adj_id with the va_id
          reconnect_root_it->second.push_back (va_adj_id);
          }
        }
      }
    }

  // initialize m_edge_pattern_index
  vector<uint> max;
  foreach (const ReconnectRoot_t::value_type& map_value, m_reconnect_root)
    {
    const uint v_id = map_value.first;
    const vector<uint>& v_adj = map_value.second;

    const uint v_degree = GraphMeasure (m_g).DegreeOf (v_id);
    const uint v_edge_pattern_max = boost::math::binomial_coefficient<double> (v_adj.size (), v_degree - m_degree);

    max.push_back (v_edge_pattern_max);
    }
  m_edge_pattern_index.Reset (max);

  // initialize m_maxNodePatterns
  m_maxNodePatterns = boost::math::binomial_coefficient<double> (m_degree-1, (m_degree-1)/2);
  }

uint 
ReconnectionPatternCounter_odd_k::GetMaxReconnectEdges () const
  {
  return m_degree - 1;
  }

uint 
ReconnectionPatternCounter_odd_k::GetMaxNodePatterns () const
  {
  return m_maxNodePatterns;
  }

vector<uint>
ReconnectionPatternCounter_odd_k::GetEdgeToReconnect (void) const
  {
  // down counter for going to a correct entry of m_reconnect_root
  uint downcounter = m_reconnect_edge_index;
  // will indicate the number of edges need to be reconnected at a node
  uint num_reconnection_edges = 0;

    // go to the node containing the edge indexed by m_reconnect_edge_index 
  ReconnectRoot_t::const_iterator va_it = m_reconnect_root.begin ();
  uint va_index = 0;
  for (; va_it != m_reconnect_root.end (); va_it++, va_index++)
    {
    uint va_id = va_it->first;
    const vector<uint>& va_adj = va_it->second;

    num_reconnection_edges = GraphMeasure (m_g).DegreeOf (va_id) - m_degree;
    if (downcounter >= num_reconnection_edges)
      {
      downcounter -= num_reconnection_edges;
      }
    else
      {
      break;
      }
    }

  // now 'it' points to a entry, where accessing to the vertex va and its neighbors (excluding the newly added vertices)
  // are available. downcounter indicates which edges to reconnect

  // get the i-th combination pattern at the node va, where r out of n edges are chosen for reconnection
  const uint n = va_it->second.size ();
  const uint r = num_reconnection_edges;
  const uint i = m_edge_pattern_index.at (va_index);
  vector<uint> perm = GetCombinationAt (n, r, i);
  const uint va_adj_index = perm[downcounter];
  uint edge_root = va_it->second.at (va_adj_index);

  // get the j-th combination pattern for the vertices {va+vb} (whose size is m), where s = m/2 vertices are choosen
  // whose neighbors are reconnected from {va+vb} to the new vertex v1  
  const uint m = m_degree - 1;
  const uint s = (m_degree - 1) / 2;
  const uint j = m_node_pattern_index;
  vector<uint> comb = GetCombinationAt (m, s, j);
  // determine if the edge is reconnected to v2
  bool isToV2 = binary_search (comb.begin (), comb.end (), va_index);
  uint edge_to;
  if (r == 1)
    {
    edge_to = (isToV2? m_v2: m_v1);
    }
  else if (r == 2)
    {
    if (downcounter == 0)
      {
      edge_to = (isToV2? m_v2: m_v1);
      }
    else
      {
      edge_to = (!isToV2? m_v2: m_v1);
      }
    }
  else
    {
    ostringstream err;
    err << "Unexpected r value - " << r << ends;
    throw exception (err.str ().c_str ());
    }

  const uint va_id = va_it->first;
  const uint ret_value[] = {edge_root, va_id, edge_to};
  return vector<uint> (ret_value, ret_value+3);
  }

bool 
ReconnectionPatternCounter_odd_k::NextEdgeToReconnect (void)
  {
  if (m_reconnect_edge_index+1 < GetMaxReconnectEdges ())
    {
    m_reconnect_edge_index++;
    return true;
    }
  else
    {
    return false;
    }
  }

bool 
ReconnectionPatternCounter_odd_k::NextReconnectionPattern (void)
  {
  m_reconnect_edge_index = 0;

  if (m_edge_pattern_index++)
    {
    m_index++;
    return true;
    }
  else
    {
    m_edge_pattern_index.Reset ();
    if (m_node_pattern_index+1 < GetMaxNodePatterns ())
      {
      m_index++;
      m_node_pattern_index++;
      return true;
      }
    else
      {
      return false;
      }
    }
  }

uint 
ReconnectionPatternCounter_odd_k::GetIndex (void) const
  {
  return m_index;
  }

void 
ReconnectionPatternCounter_odd_k::Reset (void)
  {
  m_node_pattern_index = 0;
  m_edge_pattern_index.Reset ();
  m_reconnect_edge_index = 0;
  m_index = 0;
  }

} // namespace internal
} // namespace cccn
